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The most efficient known method for solving certain computa- 
tional problems is to construct an iterated map whose fixed points 
are by design the problem's solution. Although the origins of this 
idea go back at least to Newton, the clearest expression of its log- 
ical basis is an example due to Mermin. A contemporary applica- 
tion in image recovery demonstrates the power of the method. 



1 INTRODUCTION 



Fixed points arise naturally in the study of physical phenomena, be they invariants with respect to time 
evolution (dynamical systems) or rescaling (the renormalization group). On the practical side, fixed points 
form the basis of iteration schemes specifically engineered to solve particular computational problems. One 
of the simplest applications of this idea is Mermin's solution^ 1 ) of the self-referential digit counting puzzle: 





In this sentence, 




the 


digit appears 


times; 


the 


digit 1 appears 


times; 


the 


digit 2 appears 


times; 


the 


digit 3 appears 


times; 


the 


digit 4 appears 


times; 


the 


digit 5 appears 


times; 


the 


digit 6 appears 


times; 


the 


digit 7 appears 


times; 


the 


digit 8 appears 


times; 


the 


digit 9 appears 


times. 



The object is to fill in all the blanks with decimal numerals to make the statement correct. Our instinct is to 
treat this as an exercise in logic. Mermin noted that considerably less mental energy is required by an iter- 
ative procedure, where the iterates are the vectors of ten integers which fill in the blanks: x = (xq, ... ,xg). 
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Figure 1: Left: Fourier modulus of a two dimensional object such as would be 
produced in a diffraction experiment. Right: Autocorrelation of the object obtained 
by Fourier transforming the squared Fourier modulus. 



Starting with an initial vector, say x(0) = (1, 1, 1, 1, 1, 1, 1, 1, 1, 1), successive iterates are computed by 
treating the current vector as a tentative solution and then counting the actual occurrences of the ten digits; 
thus x(l) = (1, 11, 1, 1,1, 1,1, 1,1,1). When the map encounters a fixed point, the tentative solution is 
confirmed to be an actual solution. For the choice of initial vector given above, one finds a Mermin fixed 
point already at iterate x(3) = (1, 11, 2, 1, 1, 1, 1, 1, 1, 1). As this example reveals, Mermin fixed points are 
a powerful and insufficiently explored strategy for solving a very broad range of complex problems []. 



2 IMAGE RECOVERY 

When a monochromatic plane wave is weakly and elastically scattered by an object, the resulting diffracted 
wave has an angular intensity variation derived from the Fourier modes of the object's scattering density. In 
this scenario and many others^ 2 ), access to the contrast variations within an object is available only via the 
modulus of its Fourier transform. Since phase information is usually not available, a direct inverse Fourier 
transform cannot be used to recover the object. On the other hand, relatively mundane properties of the 
object such as positivity and size, in combination with the Fourier modulus, may be sufficient to "retrieve" 
the unknown phases and recover the object. 

As an example, consider the speckle intensity pattern shown in Figure 1, formed by diffraction from an 
object in two dimensions. The data is sampled on an array measuring 448 x 448 pixels; we lack the cor- 
responding array of phases, with which a discrete Fourier transform could be used to directly produce a 

1 The reader is invited to find another Mermin fixed point of this map. There is also a 2-cycle, which solves the even harder 
problem of finding two such sentences, each of which has "this" replaced by "that." 
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Figure 2: Projections applied to the arbitrarily chosen image p WA shown on the 
left. Middle: Fourier modulus projection, 7r F (p WA ). The image shown is the mini- 
mal modification of p WA such that its Fourier transform has the modulus shown in 
Figure 1. Right. Fourier modulus projection followed by support/positivity projec- 
tion, 7T S O 7T p (p WA ). 



448 x 448 pixel image of the object. 

By the convolution theorem, the Fourier transform of the speckle gives the object's autocorrelation, shown 
on the right in Figure 1. Close examination of the autocorrelation shows its support is bounded: specifically, 
it is negligibly small outside a rectangle measuring 350 x 440 pixels. Assuming the object is positive, a well 
known theorem^ 3 ) allows one to conclude that the object's support is bounded by a rectangle measuring half 
the dimensions of the bound on the autocorrelation support, or 175 x 220 pixels. This bound on the object's 
support, we will see, is sufficient additional information to recover an image of the object. 



3 CONSTRAINTS AND PROJECTIONS 

The Fourier modulus data (Fig. 1) and the bound on the object's support are constraints that the image we 
are trying to recover must satisfy. Given an arbitrary image (for example the image on the left in Figure 2), 
the mathematical operation that restores a particular constraint, while minimizing the modification of the 
image, is called a projection^ . A natural choice of image metric is the Euclidean metric in the space of 
pixel values, since it is invariant with respect to unitary transformations into the Fourier domain. Thus the 
projection 7r P that restores the Fourier moduli consists of three operations: (1) transformation of the image 
to the Fourier domain, (2) rescaling each complex-valued pixel of the Fourier transform to the measured 
modulus (projection onto a circle), and (3) transformation back to the image domain. Another projection, 
7r s , restores positivity and the support constraint by setting pixels outside the support, and negative pixels 
within the support, equal to zero. The actions of 7r F and 7r s are illustrated in Figure 2. 

Uniqueness in image recovery requires that the number of constraints outnumber the free variables^ 5 ' 6 ) . 
Since a real- valued image of N pixels has approximately N/2 unknown phases in its Fourier transform, a 
bound on the support of the object measuring less than half the image area is normally sufficient to ensure 
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uniqueness. This condition is easily satisfied by the example introduced earlier (Fig. 1): there are 100350 
independent continuous phases in a 448 x 448 pixel Fourier transform, while the bound on the object support 
constrains a larger number of pixels, 448 x 448 — 175 x 220 = 162204, to be zero. 



4 THE DIFFERENCE MAP 



Short of having a projection that directly recovers the image by simultaneously restoring, from an arbitrary 
input, both the Fourier modulus and support/positivity constraints, one can hope to use the projections ir F 
and 7r s in an iterative fashion such that the solution can be extracted from an appropriate Mermin fixed point. 
One such approach is the difference map( 6 > : 

P^D(p) = p + /3 [lT s O f F (p) -7T F o/ s (p)] . (1) 

The action of D is to add to the current iterate p the difference of projections (composed with two additional 
maps) scaled by a parameter f3. To see how a Mermin fixed point of D, p* , provides the solution, p so \, we 
observe that D{p*) = p* implies the difference of projections vanishes. In other words, the same image, 
now identified as p so i, was produced by each of the two projections and therefore satisfies both sets of 
constraints: 

% ° /f(p*) = Psol = vr P o f s (p*) . (2) 

The maps f F and / s have so far not been specified but must be chosen with care in order to make the Mermin 
fixed point attractive. Reference 6 makes the choice 

/f = (1 + 7f)vt p - 7 F 
fs = (1 + 7s)^s - 7s , 

and finds 7 F = 7 S = — as optimal parameter values 

The difference map is superior to the naive alternating projection map A = ir s o 7r P because of stagnation 
caused by fixed points of A which do not satisfy the Fourier modulus constraint (and have no simple re- 
lationship to the solution which does, in contrast to eq. 2). To overcome stagnation, Fienup^ 7 ) introduced 
the hybrid input-output map which, interestingly, is obtained as a special case of the difference map for the 
parameter values 7 F = (3 , 7 S = — 1. Although the hybrid input-output map has been the main tool for 
phase retrieval for nearly twenty years, its fixed point properties in the geometrical setting of projections has 
come to light only recently^ 6 ' 8 ) . 



5 A PHASE RETRIEVAL EXAMPLE 



When implemented with the fast Fourier transform, the difference map can be computed in a time that grows 
only quasi-linearly with the number of pixels in the image. Although there is as yet no comprehensive 
theory of the number of iterations required to reach a Mermin fixed point, numerical experiments indicate 
that progress toward the solution, measured by the norm of the difference, 

e = Iks f F (p) ~ ttf o f s (p)\\ , (4) 
is systematic though not strictly monotone. 

2 When /p and / s are identity maps (7 F = 7s = — 1) the Mermin fixed point is found to be repulsive^ 6 ' . 
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Figure 3: Plot of log e (eq. 4) for the first 500 iterations of the difference map using 
the diffraction data of Figure 1, the support/positivity constraint as illustrated in 
Figure 2, and (3 = 1.15. 



Figure 3 shows the logarithm of the error estimate e as a function of iteration for the example introduced in 
Section 2. A rectangular object support constraint was imposed, with dimensions (175 x 220) determined 
from the autocorrelation; p WA (Fig. 2) was the initial image and (3 = 1.15 was used for all 500 iterations. 
The evolution of the iterates p(i) is shown in Figure 4. The stationarity of the last iterate, together with the 
smallness of the corresponding error (Fig. 3), leaves no doubt that a Mermin fixed point has indeed been 
found. According to equation 2, the object is recovered by applying the map 7r s o f F (or ir F o / s ) to the 
final iterate; the result is shown in Figure 5. Finally, although the fixed point p(oo) (last image in Fig. 4) 
is not unique and dependent on the initial image p(0), the uniqueness of the recovered object (Fig. 5) has 
convincingly been demonstrated over a period spanning 67 years. 



6 CONCLUSIONS 

Mermin's example^ 1 ) is the inspiration for iterative solutions to problems considerably beyond Newton's 
square root and descendants. Phase retrieval belongs to the class of feasibility problems, usually posed 
in the context of linear programming with convex constraints® . Since the Fourier modulus constraint is 
nonconvex, standard algorithms either do not apply or have no guarantee of convergence. The iterative 
difference map algorithm® , also without a bound on the number of iterations, is currently the most efficient 
method for solving a large class of phase retrieval problems. This includes a highly simplified version of 
phase retrieval: the problem of recovering a fixed length binary sequence from its cyclic autocorrelation® . 
The latter appears to be comparable in difficulty, and has a mathematical kinship to, the problem of factoring 
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Figure 4: Iterates p(l), p(5), p(50) and p(500) of the difference map; the initial 
image was p(0) = p WA (Fig. 2). 
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Figure 5: The recovered image (ir s o f F applied to the last image in Figure 4) and 
detail, after inversion. 

integers. |^| 
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3 One is interested in factoring elements in the ring of cyclotomic integers, with each of the two factors known a priori to have 
binary coefficients and related as algebraic conjugates. 
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